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Abstract 

The linear analysis of the instability due to vertical shear in the dust layer of the solar 
nebula is performed. The following assumptions are adopted throughout this paper: (1) The 
self-gravity of the dust layer is neglected. (2) One fluid model is adopted, where the dust 
aggregates have the same velocity with the gas due to strong coupling by the drag force. (3) 
The gas is incompressible. 

The calculations with both the Coriolis and the tidal forces show that the tidal force has 
a stabilizing effect. The tidal force causes the radial shear in the disk. This radial shear 
changes the wave number of the mode which is at first unstable, and the mode is eventually 
stabilized. Thus the behavior of the mode is divided into two stages: (1) the first growth 
of the unstable mode which is similar to the results without the tidal force, and (2) the 
subsequent stabilization due to an increase of the wave number by the radial shear. If the 
midplanc dust /gas density ratio is smaller than 2, the stabilization occurs before the unstable 
mode grows largely. On the other hand, the mode grows faster by one hundred orders of 
magnitude, if this ratio is larger than 20. 

Because the critical density of the gravitational instability is a few hundreds times as large 
as the gas density, the hydrodynamic instability investigated in this paper grows largely be- 
fore the onset of the gravitational instability. It is expected that the hydrodynamic instability 
develops turbulence in the dust layer and the dust aggregates are stirred up to prevent from 
settling further. The formation of planetesimals through the gravitational instabilities is 
difficult to occur as long as the dust/gas surface density ratio is equal to that for the solar 
abundance. 

On the other hand, the shear instability is suppressed and the planetesimal formation 
through the gravitational instability may occur, if dust/gas surface density ratio is hundreds 
times as large as that for the solar abundance. 
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I. INTRODUCTION 

Planetesimal formation in the solar nebula is one of unresolved issues in the theory of planetary 
formation. In the solar nebula, micron-sized dust particles which float uniformly at first are 
considered to stick each other by non-gravitational forces (e.g., the van der Waals force) to form 
dust aggregates. If the solar nebula is laminar, the dust aggregates settle toward the midplane, 
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and simultaneously they stick each other by collisions due to differential settling velocities. 
Thus, a thin dust layer is formed around the midplane. The dust density around the midplane 
increases gradually due to the dust settling. When the dust self-gravity exceeds the tidal force 
of the central star, that is, the density exceeds the critical density p c , the dust layer becomes 
gravitationally unstable and fragments into a lot of km-sized planetesimals (Safronov 1969, 
Goldreich and Ward 1973, Coradini et al. 1981, Sekiya 1983). Formerly the above mentioned 
gravitational instability model was considered to be the most promising process of planetesimal 
formation. 

However, if the solar nebula is turbulent, dust aggregates are stirred up, and it may be 
difficult for the dust density to exceed the critical density of the gravitational instability. Several 
causes of turbulence in the solar nebula have been proposed: the thermal convection (Lin and 
Papaloizou 1980, Cabot et al. 1987a, b, Ruden et al. 1988), the hydrodynamic instability due to 
radial shear in the disk (Papaloizou and Pringle 1984, Goldreich and Narayan 1985, Narayan et 
al. 1987, Nakagawa and Sekiya 1992, Sekiya et al. 1993), and the magneto-rotational instability 
(Balbus and Hawley 1991), and so on. This work concentrates on investigating the instability 
due to vertical shear in the dust layer, which arises by the dust settling itself as explained in the 
following. If the drag force did not work, dust aggregates would revolve with the Kepler velocity 
balancing the centrifugal force and the gravitational force of a central star. On the other hand, 
a gas fluid would revolve slightly slower than the Kepler velocity because only gas is partially 
supported by a gas pressure gradient. Thus, when gas and dust drag each other, their velocities 
depend on the dust to gas mass ratio and the frictional time. In particular, when the frictional 
time is much shorter than a Kepler period, namely that dust aggregate is small and/or fluffy, 
a dust aggregate and gas fluid move approximately with a same velocity. Then, the revolution 
velocity v in the co-ordinate system which moves with the local Kepler velocity vk is given by 

v = --^rjv K - (1) 

Here p is the total fluid density defined by 

P = Pg + Pd, (2) 

where p g is the gas density and p^ is the dust density. The non-dimensional parameter r] 
represents the effect of the radial pressure gradient: 

r dP 

V ~ 2p g vj c dr' {6) 

where r is the distance from the rotation axis of the disk and P is the gas pressure (Adachi et al. 
1976, Weidenschilling 1977, Nakagawa et al. 1986, Sekiya 1998). This vertical shear may cause 
the shear instability, which may develop turbulence in the dust layer; thus dust aggregates may 
be stirred up and dust density may not exceed the critical density of the gravitational instability 
(Weidenschilling 1980). 

Cuzzi et al. (1993) have performed numerical simulations of two fluids, gas and dust, in 
order to examine whether the above mentioned turbulence is strong enough to prevent the 
gravitational instability of the dust layer. They assumed single-sized compact aggregates (10cm 
or 60cm). They obtained the quasi-equilibrium density distributions balancing the turbulent 
diffusion and the settling of aggregates. They showed that the density does not reach the 
critical density of the gravitational instability. Champney et al. (1995) extended this model to 
multi- fluids where different sized dust aggregates existed in the gas. Dobrovolskis et al. (1999) 
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made calculations including the dissipation due to the friction between dust aggregates and the 
gas. The results of these papers are almost same as Cuzzi et al. (1993). The parameters used in 
the turbulence model of these papers were estimated from the laboratory experiments without 
the density stratification and the effects of the Coriolis and the tidal forces. 

Sekiya (1998) obtained the quasi-equilibrium distributions of dust density, including the 
effect of the density stratification, but neglecting the effects of the Coriolis and the tidal forces. 
The vertical shear rate dv/dz increases as dust aggregates settle toward the midplane, where z 
is the distance from the midplane. Eventually, the shear rate exceeds the critical value of the 
shear instability, i.e., the Richardson number becomes less than 0.25. The Richardson number 
is defined by 

_ g z dp/dz 

P {dv/dzf y ' 

where g z is the vertical gravitational acceleration. Then turbulence may be developed. If dust 
aggregates are small enough, very weak turbulence can stir up the dust aggregates and prevent 
them from settling further. Thus, the dust density distribution regulates itself so that the shear 
rate is marginally above the critical value J c = 0.25. Sekiya (1998) calculated analytically this 
density distribution with the critical shear rate. He has concluded that the dust density hardly 
exceeds the critical density of the gravitational instability, as long as the dust to gas mass ratio 
inferred from the solar abundance is assumed. He has shown that if the dust to gas mass ratio 
increases one order of magnitude by gas escape or dust concentration, the dust density at the 
midplane exceeds the critical value of the shear instability, and planetesimal can be formed 
through gravitational fragmentation of the dust layer. 

Youdin and Shu (2002) also found the critical ratio of the surface density ratio of dust and 
gas using the density distribution with the critical shear rate. In addition, they suggested a dust 
concentration mechanism through radial drift of chondrule-sized dust by gas drag in the nebula. 
This mechanism enables the density ratio of dust and gas to exceed the critical value in < few 
xlO 6 year. 

The derivation of dust density distribution by Sekiya (1998), and Youdin and Shu (2002) 
was somewhat intuitive. More elaborate calculations of the shear instability are needed in order 
to conclude that the gravitational instability model of planetesimal formation is denied. For 
the first step, Sekiya and Ishitsu (2000, 2001) (Hereafter Papers I and II) investigated the shear 
instability of the dust layer using the linear analysis, neglecting the effects of the Coriolis and 
the tidal forces. 

In Paper I, we assumed that the unperturbed density had the constant Richardson number 
density distributions by Sekiya (1998). The results show: (A) The flow is stable for the Richard- 
son number J £ 0.22. (B) The growth time of the shear instability is much longer than the 
Kepler period, as long as the Richardson number J <; 0.1. On the other hand, the Coriolis and 
the tidal forces would affect the flow in time scale on the order of the Kepler period. Thus the 
neglect of these forces is not good for the constant Richardson number density distribution with 
J £0.1. 

In Paper II, the linear stability analysis was performed using the following density distribu- 
tion which will be called the hybrid density distribution (an inner constant density distribution, 
plus outer sinusoidal transition regions): 



Pdo(z) 



Pdo(0) for \z\ <z d - 2h d , 
Pdo(0){l - sin[7r(z - z d + h d )/2h d }}/2 for z d - 2h d < \z\ < z d , (5) 

for z d < \z\, 
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where Zd the half-thickness of the dust layer, and hd the half-thickness of the transition zones, 
where the dust density varies from Pdo(0) to sinusoidally. Here the half-thickness of the dust 
layer is given by 

and the surface density of the dust is given by 

+°° _ / 7.1/ d (r/AU)- L5 g cm 2 for r < 2.8AU, 



oo 



30/ d (r/AU)- i - 5 g cm 2 for r > 2.8AU, 



where pd is the dust density defined by the total dust mass floating in a unit volume, and fd is 
a parameter (fd = 1 for the Hayashi model). We used Hayashi's solar nebula model (Hayashi 
1981, Hayashi et al, 1985) at 1AU as the dust surface density in the most of calculations 
except for Figs. |, ^(| and [H]. 

In the following, we explain the reason why the hybrid density distribution is used. Dust 
particles which are distributed uniformly at first stick together to form dust aggregates. In a 
laminar nebula, the settling velocity Vd z of a dust aggregate is given by 

v dz = -T f Q, 2 K z, (8) 

where Tt is the frictional time of the dust aggregate. Dust aggregates grow faster in regions with 
larger |z|, since the principal relative velocity of dust aggregates is induced by difference of set- 
tling velocities of dust aggregates with different frictional times (Weidenschilling 1980, Nakagawa 
et al. 1981). As dust aggregates grow, their settling velocities increase if the dust aggregates are 
compact. Thus, dust aggregates accumulate in a certain region with an intermediate value of \z\ 
(see 1000 yrs and 1300 yrs density distribution in Fig. 2 of Nakagawa et al. (1981)). This state 
is unstable for the Rayleigh- Taylor instability, and the dust density distribution is considered 
to be adjusted as to be constant in the dust layer (Watanabe and Yamada 2000). Thus, the 
constant distribution in the dust layer as given by Eq. (|5|) for \z\ < Zd — 2hd is considered to 
be a natural consequence of the dust settling in a laminar disk. For simplicity, the sinusoidal 
density distribution is assumed in the transition region (zd — 2hd < \z\ < Zd). 

According to results of Paper II, if pd(0)/p g ~ 1, the growth rate of the instability is on the 
order of the Keplerian angular frequency. On the other hand, if pd(0)/p g >> 1, the growth rate 
is much larger than the Keplerian angular frequency. Thus, we have expected that the Coriolis 
and tidal forces might not have an important effect as long as pd(0)/p g >> 1. 

Ishitsu and Sekiya (2002) (Hereafter Paper III), carried out calculations taking the Coriolis 
force into account, but neglecting the tidal force. The results showed that the Coriolis force 
had little effects. However, we found a new type of resonances which resembles the Lindbland 
resonances if the growth rate is similar to or smaller than the Keplerian angular frequency. The 
energy source of the instability is this resonance and different from that of the shear instability. 
These results will be reviewed in Section 4. 

In this paper, we solve linearized perturbation equations with the tidal and the Coriolis 
forces under the following assumptions: (1) The self-gravity is neglected. (2) A mixture of gas 
and dust is treated as one fluid, which is a good approximation in the case where dust aggregate 
sizes are small (^ 1cm). (3) The gas is incompressible since the dust layer which is much less 
than vertical scale height of the gas disk H is treated. (4) The effects of the radial density 
and the pressure gradient of the unperturbed state are only incorporated in the unperturbed 
rotation velocity distribution Vq(z). (5) Local Cartesian coordinates (x,y,z) are used and we 
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neglect the curvature of a circle with constant values of r and z. (6) The hybrid dust density 
distribution used in Papers II and III is adopted as the unperturbed dust density distribution. 

When the dust layer with the hybrid dust density distribution becomes unstable for a linear 
perturbation, the linear mode grows. Eventually, nonlinear effects become important, and a 
quasi-equilibrium state between settling and turbulent diffusion of dust aggregates would be 
realized. This paper treats the first linear growth of the instability. Of course, it is very 
important to obtain quasi-equilibrium state with taking account of the Coriolis and the tidal 
forces, which is left in future works. Note that J =const distribution obtained Sekiya (1998) is 
not necessarily realize, since it neglects the effects of the Coriolis and the tidal forces. 

In Section II, the basic equations for the linear analysis are derived. In Section III, a nu- 
merical method is described. In Section IV, solutions with the Coriolis force but without the 
tidal force (Paper III) are reviewed. In Section V, solutions with both the Coriolis and the tidal 
forces are given. In Section VI, conclusions obtained in this paper are summarized. 



II. FORMULATION 

We use the local Cartesian coordinate system around a radius R from the central star rotating 
with the Kepler angular frequency Qk{R)'- 

x = r-R, (9) 

y = R[<f> - n K (R)t], (10) 

z = z, (11) 

where the curvature of r = constant circle is neglected since 0(d/dr) » 0(l/r) and (d/d(f>)/r m 
(d/dy) for r >> \x\, \y\. 

In the rotational system that revolves around a central star with angular frequency Qk(R), 
the Coriolis and the centrifugal forces emerge. The tidal force, resultant force of gravity of central 
star and centrifugal force, is given in the local approximation regime by 3x0,^ (R). Hereafter, 
£Ik(R) is denoted by Qk for simplicity. Note that Qk is constant. The hydrodynamic equations 
are given for the local Cartesian co-ordinates which move azimuthally with the Kepler velocity 
vk and rotate with the Kepler angular velocity Q.k by 

du dv dw _ g 
dx dy dz 

dp ^dp ^dp ^dp _ q 
dt dx dy dz 
du du du du ldP 

at + u ai + "a» + "& = ~, * + + 2!i "'"' (14) 

dv dv dv dv 1 dP „ , . 

— + u—+v— + w— = --—-2n K u, 15 
dt dx dy dz p dy 

dw dw dw dw 1 dP ^ 2 Qg\ 
dt dx dy dz p dz K ' 

Note that only the case T = 1 is realistic; on the other hand, the case T = corresponds to 
the model used in Paper III, where the tidal force is neglected. In this paper, calculations are 
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performed for < T < 1, in order to elucidate the effect of the tidal force on the instability in 
the dust layer. If u = P = 0, Eq. ( |l4j ) reads 

v = ~Tn K x. (17) 

This expresses the circular Kepler motion in the local coordinate system (for T = 1). In order to 

eliminate the Keplerian part of the velocity, we introduce the velocity relative to the Keplerian 
motion. 

3 

v = v + -TQ, K x. (18) 

From Eqs. flT|) to (||), we have 

du dv dw 
dx dy dz 

| + „g + (s _|™ xl) | +M |.„, (20) 

du du . .du du ldP „ , . 

Tt +U d-x +{ "- 2™ KX) dy- +W Tz= -pte + 2QK *> (21) 

dv dv , .dv dv ldP 3 , . 

at + + (c " 2™*% + »& " w " (2 " 2 T)f!,f (22> 

dw dw , 3 m ^ . dw dw ldP ^ 9 , . 

«■ + + (c - 2™^ % + » & = "p & - ^ < 23 > 

In order to carry out linear calculations, we assume that an unperturbed state is steady and 
uniform in x and y directions: 

d d d 
dt dx dy 

We also assume that the unperturbed velocity is parallel to y-axis (i.e. azimuthal direction) 

uo = wq = 0. (25) 

From Eqs. (§l|), (H) and (||), we have 

1 dP 

-—^ = 2n K v , (26) 
Po dx 

and 

1 r)Pr. _ 

-ffe (27) 



Po <9z 
respectively. 

The unperturbed azimuthal velocity is calculated from a given dust density distribution p^o 
and a given value of r\ by using Eqs. (|l|) and (0): 

w = (28) 
Po 

where values of 77, vr: are assumed to be constant, and /?o depends only on z and is indepen- 
dent of x in the local approximation with r ~ i?. The midplane gas density /? s = p ff (0) is given 
by 

Pg (0) = ^E g /H, (29) 
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where S g is surface gas density: 

S 9 = 1700/ 9 (r/AU)- L5 g cm 2 , (30) 

and f g is a parameter (f g = 1 for the Hayashi model). The radial pressure gradient dPo/dx is 
then given by Eq. (p6|). 

Linearizing Eqs. (|19|) to ( |2^ ) and using Eqs. (|^) and (p7|), we have 

dui dv\ dw\ 

dx dy dz ' 

dpi 3 dpi dpo n , qo ^ 

-dT + ^-- 2 Tn KX )— + — Wl = 0, (32) 

+ _ ItQkx)^- = --^ + 2Q K % 1 + 20^, (33) 
ot 2 oy po ox po 

dvi , 3 m ^ . dvi dvn 1 dPt , 3 m , „ . , N 

— i + («o - -mu-a;)-^ + -^wa = — 1 - (2 - -T 34 

ot 2 ay dz p oy 2 

<9wi 3 ,<9wi 1 9Pi ftf^ . , 

— — + (u - -Tfi^x — — = pi. (35) 

at 2 <9y po po 

The normal mode analysis (i.e. the Fourier transform) for x cannot be done in this coordi- 
nates system because the coefficients of equations depend on x. Thus we transform coordinate 
y into the shearing coordinate y' with velocity — ^TU^x in order that the normal mode analysis 
for x can be done (see e.g., Goldreich and Lynden-Bell 1965, Goldreich and Tremaine 1978, Ryu 
and Goodman 1992): 

y' = y + ^Tn K xt. (36) 
We assume that perturbed quantities are written 

fi(x,y,z,t) = f 1 (z,t)exp[i(k y y' + k x x)]. (37) 



Substituting this expression into Eqs. (31)-(^), the perturbation equations are rewritten 

ik' x u\ + ik y v\ + -g— !- = 0, (38) 

— + ikyVopi + — w\ = 0, (39) 

+ ikyvom = -ik' x —Pi + 2Q K —pi + m K vi, (40) 

Ot po po 

^ + ikyVQVi = -iky—^P\ - ^-Wl ~ ( 2 ~ ( 41 ) 

+ zfcjpvowi = q pi, (42) 



where, 



<9i y po po 

4 = + ^Tk y n K t. (43) 

Here, we cannot use the Fourier transform method with respect to t because the coefficients of 
these equations depend on t. The linearized equations must be integrated numerically by using 
the finite difference method in order to see the growth of an instability. 
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Boundary conditions for z-direction are given as follows. Only odd solutions for w\ {i.e. 
w\ = at z = 0) are considered since even ones are always stable according to our calculations. 
The continuity of pressure at the boundary between the dust and gas layers were applied in the 
analytical method of Papers I to III, but it is difficult to apply this condition to our numerical 
method. Thus, we solve the perturbation equation numerically within region [0, zq] where the 
solid-wall condition is applied at the boundary for simplicity. The value of z$ should be large 
enough in order for an eigenfunction to decay sufficiently at a boundary zq. It was confirmed 
that numerical solutions for T = with this condition agree well with solutions of the Paper 
III. Thus, the boundary conditions are given by 

w\ = at z = and zq. (44) 

From Eqs. ( p^ ) and (|4"4|), we have 

<9Pi 



dz 

From Eqs. ( |39]) and d44|), we have 

dp 



+ Vt K zp\ = at z = and zq. (45) 



— ^ + ikyVQpi = at z = and Zq. (46) 



If we give an initial condition with p\ = at z = and zq, Eq. ( |4q ) implies that p\ = for t > 0. 
Hereafter we take only such initial conditions for simplicity. Accordingly boundary conditions 
for Pi and p\ are given by 

dP 

— — = at z = and zq, (47) 
dz 

pi = at z = and Zq. (48) 



III. NUMERICAL METHOD 

We adopt MAC method (Harlow and Welch 1965), i.e. pressure is determined by the con- 
dition that the equation of continuity is satisfied in the next step, and other variables ui,vi,wi 
and p\ are calculated using this value of pressure. 

We define the divergence of a perturbed velocity by 



D = ik' x m + ikyVi + — ^. (49) 



dz 

Multiplying Eq. (|4^) by ik' x and Eq. ( |4l| ) by ik y , and taking partial derivative Eq. ( |42| ) with 
respect to z, and adding, we have 



9D _ .-L^n o;u d *°.... 1 7/2 .2, ^ 1 d. i ldpo9fl 



at " lhydoD ~ 2iky lz~ m ~ 7 {-** - k y + dz^) Pl + pidz dz 

+2in K k' x —p 1 + 2in K k' x vi - i (2 - 3T) Q K k yUl + Q 2 K — - (-—pi ) ■ (50) 
Po dz V Po 



We use the first order approximation: 



dD _ L> n+1 - D n 
~dt ^ A* ' 



(51) 
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We require that the equation of continuity (|38|) is satisfied in the next step, i.e. D 
Thus, we get the Poisson equation for Pi: 

( h m2 k 2,d 2 1 dp d\ 

' n.-j. - \ Tin o.-?. _ _ n | ,.//),/) 



n+1 



po \ ( -r- - 2zfcj,u ) - 2ik y —^w^ + 2in K k'™v^ - % (2 - 3T) n^u? 



+2 i n K k' x ^p n 1 + Q^-f- (_ A p n ) l . (52 ) 
Po V Po /J 

From this equation, we can get Pi at the time step n. Next, we calculate pi,u\,vi and w\ at 
the time step n+1 from Eqs. (|39|)-(^): 

dpo„„„ 
dz 

u™ +1 = v% + At l-ikyVQU^ - ik'^—P? + 20 K %? + 2^<1 , (54) 
I Po Po J 

v1+ l = vl + At j-^ < - - »Jfe„— i? " ( 2 " } » ( 55 ) 

= < + Atl-ikyvow? - - ^p?) • (56) 

t Po <9z Po J 

The above method is the first order accuracy with respect to At. In order to improve them 
to the second order accuracy, we use the following strategy. We replace perturbed quantities 
fi = (Pi ) u i > v i 3 ^i ) on right hand side of Eq. ( |52"D with 

/f 1 = + /i)/2- (57) 

Again, we solve Eq. ( |52"| ) replacing n with n + 1/2. More exact values of + , + , v™ + and 
are obtained by replacing these quantities at n in braces of the right hands of Eqs. ( |53| ) 
to ( p6[ ) with ones at n + 1/2. We adopt one dimension staggered grid, where w\ are estimated 
at grids and Pi, pi,ui,v\ at midpoints of adjacent grids. In Eqs. (|52|), (|5^) and fl55|) , u;i is 
calculated by taking the mean values at the adjent meshes. So as p\ in Eq. (|D 



=p n 1 +At{ -ikyvofi ~ -fV , (53) 



Numerical parameters (Atftx, Zo/rjr, Grid number for z-direction N z ) and model parameters 
(pd(Q)/p g , hd/zj, fg, fd, k x r]r, log(kyT] 2 r 2 ), T) in this work are listed in Table I. The minimum 
values of the Richardson number J m in for a region [0, zq[ are also listed as an indicator of shear 
strength. 

Initial conditions were set as follows. Some Fourier components of lower orders were selected 
as to satisfy boundary conditions of tii, w\ and p\. Each Fourier coefficient is given by a random 
number. Velocity v\ was determined from u\ and w\ by using the equation of continuity (^). 

IV. SHEAR INSTABILITY WITHOUT THE TIDAL FORCE 

Here, solutions without the tidal force, that is, the case T = is obtained in order to compare 
the results with the tidal force (detailed description is given in Paper III). Then, we can carry 
out the normal mode analysis for time: 

f 1 (z,t) = f 1 (z)exp(-iwt), (58) 
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where u is the angular frequency of perturbed quantities, and obtain eigenvalues oj by solving 
the differential equations to satisfy with boundary conditions. The growth rate of the instability 
is ui = 9(cj), where denotes the imaginary part. 

Figure |l] shows the growth rate oji as a function of radial and azimuthal wave number k x and 
k y in the case where Pdo(0)/p g = 0.5, and hd/zd = 0.2. The growth rate uji has a finite positive 
value for k x < k xc , and uj = for k x > k xc , where k xc is the critical radial wave number. This 
is very important for the stability with the tidal force as will be written in the next section. 

Figure || shows the growth rate of the instability with the most unstable wave number 
(hereafter called "the peak growth rate") as a function of Pd(fy/Pg with h^/zd = 0.2. Note the 
slope of curve in Fig. || changes at pd(0)/p g = 0.23 and 0.39. This is because one mode which 
has the peak growth rate at smaller dust density at the midplane moves to another mode at the 
two dust densities at the midplane. The increase of the dust density corresponds to the dust 
settling to the midplane in a quiescent nebula. Thus, in the case without the tidal force, the 
growth rate of the shear instability increases as the settling proceeds. The peak growth rate also 
increases when the half-thickness of the transition zone of the dust layer hd decreases (Fig. ||). 

Figure ||| shows the peak growth rate as a function of fd (see Eq. (^)). As fd increases, the 
peak growth rate decreases. When fd increases with a constant dust density at midplane, the 
half-thickness of transitional zone of the dust layer increases. This implies that the shear rate 
dvo/dz decreases. Thus, the shear instability is depressed if the dust surface density increases 
due to some mechanism. 

Figure |5| shows the peak growth rate as a function of f g (see Eq. (j30|)). In contrast with fd, 
as f g decreases, the peak growth rate increases. Figure |] is identical to Figure || when 1/ f g is 
taken as fd- Thus, only will be taken as a numerical paramter in the following. 



V. SHEAR INSTABILITY WITH THE TIDAL FORCE 

We performed a test run for the case T = in order to check our code, since we have already 
gotten eigenfuctions of growing modes in the previous section. Figure |6| shows time evolutions of 
the radial, azimuthal and vertical components of the mean kinetic energy for log(fc^ 2 r 2 ) = 2.46, 
T = 0, pdo(0)/p g = 0.5 and hd/zd = 0.2. Note that important values are not the absolute values 
of components of the perturbed kinetic energy but ratios and the growth rate of those because we 
treat linearized quantities. If a linearized variable grows proportional to exp((j/t), its absolute 
value squared varies as f\f* ex exp(2u>/t), where the superscript * denotes the complex conjugate 
value. Thus, we calculate the growth rate by 



u>r = — — In 
2At 



(< uxuX > + < v\v{ > + < w\w\ >)\t+M 



(59) 



(< U\U\ > + < V\V{ > + < wxw* x >)\ t 

The time evolution of the growth rate is shown in Fig. ^. We find that the growth rate 
asymptotically approaches the analytical solution (expressed by the dotted line). In addition, 
we have checked that numerically calculated function w\(z, t) for large values of t is almost same 
as the eigenf unction. Thus we consider our code works correctly. 

Figure || shows the time evolutions of radial, azimuthal and vertical components of mean 
kinetic energy with same parameters as Fig. || except that T = 1. The solid curve in Fig. || 
shows the temporal growth rate with log(/c 2 7/ 2 r 2 ) = 2.46. When T = 1, the growth rate rises 
at first and then begins to vibrate, being zero in average. Thus the tidal force has a stabilizing 
effect. 

The dotted curve in Fig. |9] shows the temporal growth rates with log(fc 2 r/ 2 r 2 ) = 2.00. The 
temporal growth rate is rather insensitive to the wave number. The solution with a higher wave 
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number has a tendency to transit to simple vibration earlier. Initial conditions have influence 
on growth of the energy at the beginning. However, subsequently, time evolution of the growth 
rate and energy become independent of initial conditions. 




Figure 10 shows the growth rate with different values of the tidal parameter T. Note that 
T = 1 is only the realistic case, and we solved cases with other values of T for comparison in 
order to elucidate the effect of the tidal force. As T increases, the both amplitude and period of 
oscillation of ujj increase. We have checked that these oscillation frequencies are almost equal 
to the epicyclic frequency, i.e. the oscillation frequency due to the Coriolis and the tidal forces 
(see e.g. Eq. (8.9) of Shu 1992) 

(60) 

and are different from the tidal shear rate 3T$1k/2- Here we consider the reason why the 
instability is stabilized even if T is as small as 0.1. As found from Eq. (f43|), k' x increases as 
time t elapses. The analysis for T = indicates that, as k x increases, the growth rate of the 
instability decreases and eventually the instability disappears (see Fig. |]). Consequently for 
T / 0, even if an unstable mode had grown at the beginning, it would be stabilized as time 
elapses. A time for stabilization t s is estimated as follows. Given the critical wave member k xc 
above which tuj = for the case of T = 0, we have 

t s = {kxc - k x )j (^k y ft K T\ , (61) 



from Eq. ([43|). It is understood from Eq. (31) that the time scale to stabilize the instability 
decreases as wave number k y increases. The oscillations of kinetic energy and the growth rate 
are interpreted as mutual interference of non- growing waves with ivj = 0. We have found that 
the oscillation period hardly change when k x > k xc in the case of T = 0. Indeed, the oscillation 
period for Pdo(fy/Pg = 0.5, hd/zd = 0.2, T = and k x r]r = 50 becomes almost same as that 



for T = 1 as time elapses (compare Figs. |8| and 11), although the oscillation for T = 1 seems 
different from that for T = at the early phase due to relicts of initial conditions. 

Next, we consider the case where Pdo(0)/Pg = 2. The results of the previous section show 
that the peak growth rate of the instability for Pdo{fy/Pg = 2 is much larger than the case 
Pdo(0)/p g = 0.5 as long as the effect of the tidal force is not taken into account, i.e. T = 



(see Fig. g). Figures [L2] and 13 show the time evolutions of components of energy and the 
growth rate, respectively. Reflecting analytical results in the case of T = 0, there is a steep 
rise in the early stage. The growth rate decays during about a Kepler period after the early 
increase. Although the instability grows several orders of magnitude before stabilization (Fig. 
12), we cannot definitely claim whether the non- linear effect should work or not, since the initial 



amplitudes are unknown. 

Figures [14] and [0] show absolute values of eigenfunction for Pdo(fy/Pg = 2, hd/zd = 0.2, 
fd = 1, T = 1 and \og{kyT] 2 r 2 ) = 3.55 when the radial kinetic energy the smallest (tClx = 25.6) 
and the largest (t^lpc = 27.1) values, respectively. A bold solid curve denotes the Richardson 
number distribution. The vertical velocity \wi\ is omitted because it is of two orders smaller 



than |ki| and \vi\ (see Fig 12). The eigenfunctions have the largest amplitude at the height 
where the Richardson number has the smallest value. 

From Fig. |3| for T = 0, a small initial growth rate of the instability is also expected for a 



large value of hd/zd- Figures [16| and |17J show the time evolutions of components of energy and 



the growth rate for hd/zd = 0.5, log(^ry 2 r 2 ) = 2.46, T = 1 and Pdo{0)/Pg = 2, respectively. 
This growth rate rises at first and then begins to vibrate, like the case where Pdo{fy/Pg = 0.5 
and hd/zd = 0.2. 
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From Fig. |2|, the growth rate in the case T = is obviously very large when dust set- 
tling proceeds and Pdo(fy/Pg > 20. We calculated numerically the time evolutions in the 
cases \og{k 2 y rfr 2 ) = 6.18 with T = 1, p d o(0)/pg = 20 and h d /z d = 0.2 (Fig. [l|). When 
log(k 2 rj 2 r 2 ) = 6.18, we found k xc = 2.4 x 10 3 for T = 0, and then we have t s Qx = 1-3 from 



Eq. (|6l|). The estimated value of the time for stabilization is displayed by the arrow in Fig. |18|, 
which is consistent with the stabilization time of the numerical solution. During several Keple- 
rian periods, the perturbed kinetic energy increases by more than 100 orders for pdo(0)/p 9 = 20 
in this linear calculation in the case logO^V) = 6.18 (Fig. |l|). 

Our results indicate that the perturbation increases sufficiently, so that the nonlinear effect 
would be important, before the mode is stabilized. Thus, the tidal force is effective only in 
early phase of the dust settling where Pdo(®)/Pg ~ 1- The shear instability eventually grows 
to nonlinear phase when p^(0)//9 9 becomes on the order of 10 by the dust settling as long as 
fd/fg = 1 (i.e. the dust-gas surface density ratio is given by the solar elemental abundance). 

Last, we consider the case where Pdo(®)/Pg = 20 and fd = 100 to see the dependency of 
the shear instability on fd, that is, the dust surface density enhancement relative to that the 
solar abundance (see Figs. ^ and ^). The shear instability is stabilized before the perturbed 
kinetic energy grows largely in this case. According to results without the tidal force, the peak 
growth rate in this case is about a Kepler frequency and comparable to that in the case with 
Pdo(0)/p g ~ 0.5 and fd = 1 (compare Fig. |2| and . Thus, even though pdo(0)/p g increases, the 
shear instability is stabilized and planetesimals can be formed by the gravitational instability if 
dust concentrates to some regions in the solar nebula. 

As stated in §IV, decrease of f g is idntical to increase of fd- This implies that the local 
dissipation of gas also causes stabilization of the shear instability. 



VI. CONCLUSIONS AND DISCUSSION 

The linear stability analysis of the dust layer in the solar nebula is done including the effects of 
both the Coriolis and the tidal forces. The hybrid density distribution of the dust is assumed 
where the dust density is constant in the dust layer and the transition region has a sinusoidal 
density distribution. The calculations with both the Coriolis and the tidal forces show that the 
tidal force has a stabilizing effect. The tidal force causes the radial shear in the disk. This 
radial shear changes the radial wave number of the mode which is at first unstable, and the 
mode is eventually stabilized. Thus the behavior of a mode is divided into two stages: (1) the 
first growth of the unstable mode which is similar to the results without the tidal force, and (2) 
the subsequent stabilization due to an increase of the radial wave number by the radial shear. 
If Pd(0)/p g ^ 2, the stabilization occurs before the unstable mode grows largely. On the other 
hand, the mode grows more than hundreds orders of magnitude, if p d (0)/p g ^ 20. 

Since p c /Pg is much larger than 20, the linear hydrodynamic instability calculated in this 
paper grows largely before dust settling proceeds enough for the dust layer to be gravitationally 
unstable. If this instability develops turbulence in the transition region between the dust layer 
and the gas layer, the overshoot of the turbulence would widen the transition region, and the 
transition region would be stabilized, because the growth rate of the shear instability decreases 
as the half-thickness of the transition region increases. Then, the dust aggregates settle further, 
and the transition region becomes unstable again. Repeats of this cycle would eventually lead 
the whole of the dust layer to be a quasi-equilibrium state where the dust settling and turbu- 
lent diffusion balance each other. Note that this equilibrium state is not necessarily given by 
J=constant distribution derived by Sekiya (1998), because the constant J distribution is derived 
by neglecting the effects of the Coriolis and tidal forces. If this equilibrium state is realized be- 
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fore the density exceeds the critical density of the gravitational instability, planetesimals cannot 
be formed by the gravitational instability. Further nonlinear calculations on the evolution of the 
dust layer should be done in the future in order to see whether turbulence really develops and 
an equilibrium state is realized. 

We have assumed that the surface density distribution is given by the Hayashi model (i.e. 
fd = 1, f g = 1) for most parts of calculations in this paper. Our results show that the shear 
instability is stabilized by the tidal force if the ratio of the dust surface density E<j to the 
gas surface density £ s is hundreds times as large as that calculated from the solar elemental 
abundance, even if Pdo(0) ~ p c - The shear-induced instability does not grow largely when the 
dust density reaches the critical density of the gravitational instability in this linear calculation 
if dust concentrates to some locations, or gas dissipates from the nebula. For example, dust 
concentration due to radial migration induced by the gas drag in the laminar disk could increase 
the dust to gas mass ratio (Youdin and Shu 2002). Then the planetesimal formation due to 
gravitational instability may occur. In this case, the effect of the self-gravity, e.g. midplane cusp 
of the vertical density distribution (Sekiya 1998, Youdin and Shu 2002) should be taken into 
account. 

Lastly, we discuss the dust concentration by eddies in a protoplanetary disk. Cuzzi et al. 
(2001) assumed a turbulent disk and showed that the dust concentrates to preferential positions 
between eddies, even if the dust aggregates are as small as 1mm. In their analysis, the effects of 
the solar gravity are not included. If dust concentrated to some positions, they revolve around 
the central star with the Keplerian velocity. On the other hand, the residual gas rich parts 
revolve with smaller mean velocity due to the radial pressure gradient. Thus, the dust aggregates 
may be removed from the preferential positions between eddies and the concentration may be 
suppressed. There are also models of dust concentration by large eddies. Klahr and Henning 
(1997) assumed large scale eddies whose rotation axis are parallel to the midplane. They showed 
that dust concentration away from the midplane. On the other hand, another "vortex capture" 
model (Barge and Sommeria 1995, Tanga et al. 1996, Godon and Livio 1999, 2000, Chavanis 
2000, de la Fuente Marcos and Barge 2001) assumes eddies whose rotational axis is parallel to the 
rotation axis of the protoplanetary disk. These eddies preferentially concentrate dust aggregates 
which have stopping time due to gas drag on the order of the Keplerian period. These models 
give an additional possibility to reach the gravitational instability. However, it is not clear what 
kind of eddies are really probable in the protoplanetary disks. Full hydrodynamic simulation 
including interactions between gas and dust and the solar gravity should be performed in future 
to elucidate whether the dust concentration by eddies really occurs in the protoplanetary disks. 
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Figure 1: The growth rate toj of the mode as a function of radial and azimuthal wave number 
k x and k y in the case where T = with pdo(®)/Pg = 0.5 and hd/z^ = 0.2. 
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Ao(0) 1 P g 

Figure 2: The growth rate ui of the mode with the most unstable wave number as a function of 
the ratio of dust to gas on the midplane pdo(ty/Pg, m the case where T = with hd/zd = 0.2. 
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Figure 3: The growth rate of instability oji of the mode with the most unstable wave number 
as a function of the hd/zd in the case where T = with Pdo(0)/Pg = 2. 
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Figure 4: The growth rate of instability ujj of the mode with the most unstable wave number 
as a function of the fd in the case where T = with Pdo(fy/Pg = 20 and ha/za = 0.2. 
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Figure 5: The growth rate of instability ujj of the mode with the most unstable wave number 
as a function of the f g in the case where T = with pdo(0)/Pg = 20 and hd/zd = 0.2. 



22 




Figure 6: The time evolution of radial, azimuthal and vertical parts of the perturbed kinetic en- 
ergy density, in the case where T = 0, with log(kyi] 2 r 2 ) = 2.46, pdo{0)/p 9 = 0.5 and h^/zd = 0.2. 
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Figure 7: The time evolution of the growth rate, in the case where T = 0, with 
log(ky7] 2 r 2 ) = 2.46, Pdo(ty/pg = 0.5 and h^/z^ = 0.2. The dotted line shows the analytical 
solution. 
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Figure 8: Same as Fig. || except that T = 1 and the position of a arrow shows the stabilization 
time t s flK- 
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Figure 9: The time evolution of the growth rate for log(k 2 r] 2 r 2 ) = 2.46 (i.e. the parameters ha 
same values as Fig. ||) and log(k 2 i] 2 r 2 ) = 2.00, with T = 1, Pdo{0)/p 9 = 0.5 and h^/z^ = 
The positions of solid and dotted arrows show the stabilization times t s Q,K for \og{k 2 rj 2 r 2 ) = 2. 
and 2.00, respectively. 
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Figure 10: The time evolutions of the growth rate for T = 1, 0.5 and 0.1, with 
log(ky7] 2 r 2 ) = 2.46, pdo(0)/p g = 0.5 and hj/zd = 0.2. The positions of solid, dotted, and 
dashed arrows show the stabilization times t s Q.K for T = 1, 0.5 and 0.1, respectively. 
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Figure 11: The time evolution of radial, azimuthal and vertical parts of the perturbed kinetic en- 
ergy density for T = and k x rfr = 50, with log(&^r/ 2 r 2 ) = 2.46, Pdo{0)/p g = 0.5 and hd/zd = 0.2. 
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Figure 12: The time evolution of radial, azimuthal and vertical parts of the perturbed kinetic 
energy density, in the case where h^/z^ = 0.2 with pdo{0)/p g = 2, T = 1 and \og{k 2 y r] 2 r 2 ) = 3.55. 
The position of arrow shows the stabilization time t s Q.K- 
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Figure 13: The time evolution of the growth rate, in the case where hd/zd = 0.2 with 
Pdo(0)/p g = 1, T = 1 and log^r/ 2 ?- 2 ) = 3.55. The position of arrow shows the stabilization 
time t s ^}K- 
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Figure 14: Eigenfunctions of the perturbed radial (solid curve) and azimutal velocity (dotted 
curve) at a time when < u\ > /2 reaches a minimum value and < v\ > /2 reaches a maximum 
value ( iSTft- = 25.6) in the case where pdo(0)/pg = 1, h^/z^ = 0.2, f g = 1, = 1, T = 1 and 
log(£^r/ 2 r 2 ) = 3.55. A Bold solid curve denotes the Richardson number distribution. 
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Figure 15: Same as Fig. 14 except for a time when < u\ > /2 reaches a maximum value and 
< uf > /2 reaches a minimum value ( tf2^ = 27.1). 
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Figure 16: The time evolution of radial, azimuthal and vertical parts of the perturbed kinetic 
energy density, in the case where h^/z^ = 0.5 with pdo{0)/p g = 2, T = 1 and \og(kyi] 2 r 2 ) = 2.46. 
The position of arrow shows the stabilization time t s Q,K- 
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Figure 17: The time evolution of the growth rate, in the case where hd/zd = 0.5 with 
Pdo(0)/p g = 2, T = 1 and log(kyT] 2 r 2 ) = 2.46. The position of arrow shows the stabilization 
time t s ^}K- 
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Figure 18: The time evolution of the growth rate, in the case where \og{k 2 tf'r 2 ) = 5.96 with 
T = 1, P(2o(0)/ p g = 20 and h^/zd = 0.2. The position of a arrow shows the stabilization time 
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Figure 19: The time evolution of radial, azimuthal and vertical parts of the perturbed kinetic en- 
ergy density, in the case where \og{k 2 ri 2 r 2 ) = 5.96, with T = 1, Pdo(0)/Pg = 20 and h^/zd = 0.2. 
The position of arrow shows the stabilization time t s Qx- 
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Figure 20: The time evolution of the growth rate, in the case where fd = 100 and Pdo(0)/p g = 20, 
with \og{kyT] 2 r 2 ) = 2.34 and hd/zd = 0.2. The position of arrow shows the stabilization time 
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Figure 21: The time evolution of radial, azimuthal and vertical parts of the perturbed kinetic 
energy density, in the case where fd = 100 and Pdo{fy/Pg = 20, with log(£^r/ 2 r 2 ) = 2.34 and 
hd/zd = 0.2. The position of arrow shows the stabilization time t s Q.K- 
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